/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the unobserved components model for the NKPC with stochastic trend inflation
with two structural breaks. The variability of shocks to trend inflation are constrained to be constant
across regimes.

Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k_{S_t} Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t ~iid.N(0,sig^2_v)
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_t~iid.N(0,sig^2_{u,D_t})
z_t = psi_{S_t} z_{t-1} + epsilon_t, epsilon_t|S_t ~iid.N(0,sig^2_{epsilon,S_t})

u_t and epsilon_t are correlated with correlation coefficient rho.

S_t = {0,1,2}
Prob[ S(t)=0 | S(t-1)=0] = p00
Prob[ S(t)=1 | S(t-1)=1] = p11
Prob[ S(t)=2 | S(t-1)=2] = 1

sig^2_{u,D_t} = sig^2_{u,0}(1-D_t) + sig^2_{u,1}D_t
D_t = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 

qtr= data[startsmpl:endsmpl,1];  
/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
inff1=data[startsmpl+1:endsmpl+1,8];  @one-quarter ahead inflation@

xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~xx);

@Define dummy variable series D_t to account for the Great Moderation@
GM = 130; @1984:2@
dmat = ones(GM,1)|zeros(t-GM,1); 

@================= Initialize Global Variables============@

tvp_no=4; 
start= 11;
t1=75;t2=125;

@Starting Values for PCE inflation@
prmtr_in=   4.3360830|
            3.6207972|
            1.0297822|
            0.45991641|
            0.34008988|
            0.73952899|
            1.5544761|
            0.84326407|
            -0.083434780|
            1.4732290|
            -4.3488046|
            -0.077438017|
            9.2330093|
            -0.0016211357|
            0.018851190|
            0|
            0;
 /*Robustness checks*/
 @Starting Values for CPI inflation@
 /*prmtr_in=    4.3489661 |
                3.9678287 |
                1.0310126 |
                0.45915495 |
                0.29570007 |
                1.0090050 |
                2.2200207 |
                1.1342445 |
                -0.079288579 |
                1.4791923 |
                -4.3570819 |
                0.91213739 |
                6.1668916 |
                -0.01 |
                0.018944040 |
                0.021589511 |
                0.035516920;*/  
 
  @Starting values for GDP chain-weighted inflation@
  /*prmtr_in=   4.3360830 |
                4 |
                1.0297822 |
                0.45991641 |
                0.34008988 |
                0.73952899 |
                1.5544761 |
                0.84326407 |
                -0.083434780 |
                1.4732290 |
                -4.3488046 |
                -0.077438017 |
                 9.2330093 |
                -0.0016211357 |
                 0.018851190|
                 0|
                 0.02 ;*/
  
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);
PRM_FNL=TRANS(xout);           
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl));        @Standard errors of the estimated coefficients@
 
{B_TT,PR_TT,PR_TL,stdprederr,xfor,confint,twostep} =FILTER(XOUT);
smooth3=SMOOTH(pr_tt,pr_tl,prm_fnl[1:2,1]); @Pr[S_t|Y_T], Smoothed probabilities@

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";
 
@Produces Figure A7 and A8 in the unpublished appendix@
_pdate=0;
_pltype = {5 6 3};
_plegctl={4 4 6 5};
_plegstr="Regime1\000Regime2\000Regime3";
title("Smoothed probabilities of regimes");
title("");
XY(qtr[start:t],SMOOTH3); 

@Produces Figure A9 and A10 in the unpublished appendix@
_pltype = {6 3 5 5};
_plegctl={2 4 5.5 5.3};
_plegstr="Trend Inflation\000Actual Inflation\000Lower and Upper\00090\% Confidence Bands";
title("Actual inflation and trend inflation");
title("");
xy(qtr[start:t],b_tt[.,1]~inf[start:t]~confint);

@Produces Figure A11 and A12 in the unpublished appendix@
_pltype = {6 3};
_plegctl={2 4 6 5.7};
_plegstr="Inflation Cycle\000Inflation Gap";
title("Inflation cycle and inflation gap");
title("");
XY(qtr[start:t],xfor~inf[start:t]-b_tt[.,1]);

@Produces figure A2 in the unpublished appendix@
_pltype = {6 3};
title("Actual one-quarter-ahead inflation and predicted one-quarter-ahead inflation");
title("");
_plegctl={2 4 4.3 5.8}; 
_plegstr="Actual one-quarter-ahead inflation\000Predicted one-quarter-ahead inflation";
xy(qtr[start:t],inff1[start:t]~twostep);
  
output file=UC1.out reset;
output off;
end;

@=================================================================================@

PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  
 
EXTERNAL PROC TRANS, V_PROB;
 
    PRMTR=TRANS(PRMTR1);
        
    P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
    P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
    P22=1;             @Pr[St=2/St-1=2]@
                
    siguGM0 = prmtr[3,1];
    siguGM1 = prmtr[4,1];
    sigv= prmtr[5,1]; 
    sigep = prmtr[6:8,1]; 
    rho = prmtr[9,1];
                    
    phi1 = prmtr[10,1];
    phi2 = prmtr[11,1];
    psi=prmtr[12:14,1];    
    k=prmtr[15:17,1];
    
    A=1~0;      @selection matrix@
    H = 1~1~0~0|0~0~1~0;
               
    B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
    P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
    B_LL1=inf[t1]|0|xx[1]|0; 
    P_LL1 = eye(4)*10000;       
    B_LL2=inf[t2]|0|xx[1]|0; 
    P_LL2 = eye(4)*10000;      
                
    subF = phi1~phi2| 1~0; 
    subQ= siguGM0^2~0|0~0;
    vecp = inv(eye(4)-subF.*.subf)*vec(subq);
    vv = reshape(vecp',2,2);
    p_ll0[3:4,3:4] = vv;
    p_ll1[3:4,3:4] = vv;
    p_ll2[3:4,3:4] = vv;

    p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
    p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
    p_ll2[2,2] = sigep[3,1]^2/(1-psi[3]^2);

    PROB__0 = 1; @Unconditional probability for St=0@
    PROB__1 = 0; @Unconditional probability for St=1@
    PROB__2 = 0; @Unconditional probability for St=2@
 
/*  START ITERATION */
/*===================================================================== */

    LIKV=0.0;
    J_ITER = 1;
    do until J_ITER>T;

        yt=inf[j_iter]|xx[j_iter];
        dmt = dmat[j_iter];
        sigu=siguGM0*dmt + siguGM1*(1-dmt);
 
        Q0=    sigv^2~0~0~0|
               0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
               0~rho*sigu*sigep[1,1]~sigu^2~0|
               0~0~0~0;
        Q1=    sigv^2~0~0~0|
               0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
               0~rho*sigu*sigep[2,1]~sigu^2~0|
               0~0~0~0;
        Q2=    sigv^2~0~0~0|
               0~sigep[3]^2~rho*sigu*sigep[3,1]~0|
               0~rho*sigu*sigep[3,1]~sigu^2~0|
               0~0~0~0;
        F0 =   1~0~0~0 |
               0~psi[1]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F1 =   1~0~0~0 |
               0~psi[2]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F2 =   1~0~0~0 |
               0~psi[3]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;

        FF= phi1~phi2|1~0;
        xfor = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
        ome0 =0~0~xfor| 0~0~0~0;
        xfor = k[2] *A*(inv(eye(2)-FF)); 
        ome1 =0~0~xfor| 0~0~0~0;
        xfor = k[3] *A*(inv(eye(2)-FF)); 
        ome2 =0~0~xfor| 0~0~0~0;

        /* KALMAN FILTER */
        @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
        B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
        B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
        B_TL02 = F2 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=2@
        B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
        B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
        B_TL12 = F2 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=2@
        B_TL20 = F0 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=0@
        B_TL21 = F1 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=1@
        B_TL22 = F2 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=2@
        
        @ P_tl(i,j),  S(t-1)=i, S(t)=j@
        P_TL00 = F0 * P_LL0 * F0' + Q0;
        P_TL01 = F1 * P_LL0 * F1' + Q1;
        P_TL02 = F2 * P_LL0 * F2' + Q2;
        P_TL10 = F0 * P_LL1 * F0' + Q0;
        P_TL11 = F1 * P_LL1 * F1' + Q1;
        P_TL12 = F2 * P_LL1 * F2' + Q2;
        P_TL20 = F0 * P_LL2 * F0' + Q0;
        P_TL21 = F1 * P_LL2 * F1' + Q1;
        P_TL22 = F2 * P_LL2 * F2' + Q2;
 
        @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
        F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
        F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
        F_CAST02= yt- H * B_TL02 - ome2* B_TL02;
        F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
        F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
        F_CAST12= yt- H * B_TL12 - ome2* B_TL12;
        F_CAST20= yt- H * B_TL20 - ome0* B_TL20;
        F_CAST21= yt- H * B_TL21 - ome1* B_TL21;
        F_CAST22= yt- H * B_TL22 - ome2* B_TL22;
 
        @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
        SS00= H * P_TL00 * H' ; 
        SS01= H*  P_TL01 * H' ;
        SS02= H*  P_TL02 * H' ;
        SS10= H * P_TL10 * H' ;
        SS11= H * P_TL11 * H' ;
        SS12= H * P_TL12 * H' ;
        SS20= H * P_TL20 * H' ;
        SS21= H * P_TL21 * H' ;
        SS22= H * P_TL22 * H' ;
 
        /*Updating*/
        B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
        B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
        B_TT02 = B_TL02 + (P_TL02 * H' *inv(SS02)) * F_CAST02;
        B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
        B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
        B_TT12 = B_TL12 + (P_TL12 * H' *inv(SS12)) * F_CAST12;
        B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
        B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
        B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

        P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
        P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
        P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
        P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
        P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
        P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
        P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
        P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
        P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;
       
        /* HAMILTON FILTER */   

        PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
        PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
        PR_VL02=0;
        PR_VL10=0;
        PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
        PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
        PR_VL20=0;
        PR_VL21=0;
        PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
        
        PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                                  
        @Pr[St,St-1/Yt]@
        PRO_00=PR_VL00/PR_VAL;       
        PRO_01=PR_VL01/PR_VAL;
        PRO_02=PR_VL02/PR_VAL;
        PRO_10=PR_VL10/PR_VAL;
        PRO_11=PR_VL11/PR_VAL;
        PRO_12=PR_VL12/PR_VAL;
        PRO_20=PR_VL20/PR_VAL;
        PRO_21=PR_VL21/PR_VAL;
        PRO_22=PR_VL22/PR_VAL;
        PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
        PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
        PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@

       /* COLLAPSING TERMS */
      
       if(prob__0!=0);
       B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
       endif;
       if(prob__1!=0);
       B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
       endif;
       if(prob__2!=0);
       B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
       endif;
       if(prob__0!=0);
       P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
       PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
       PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
       endif;
       if(prob__1!=0);
       P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
       PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
       PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
       endif;
       if(prob__2!=0);
       P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
       PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
       PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
       endif;
       if J_ITER < START; goto skip; endif;
       LIKV = LIKV -LN(PR_VAL);
       skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
  
RETP(LIKV);
ENDP;
 
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @constraining values of reg. coeff.@
local prm_new,d0,d1;

prm_new=prmtr1; 

prm_new[1:2,1]=EXP(prmtr1[1:2,1]) ./ (1 + EXP(prmtr1[1:2,1])); 
prm_new[3:8,1] = sqrt(prmtr1[3:8,1]^2);
prm_new[9] =prmtr1[9,1]./(1+abs(prmtr1[9,1]));
prm_new[12:14,1] =prmtr1[12:14,1]./(1+abs(prmtr1[12:14,1]));
d0 =prmtr1[10,1]./(1+abs(prmtr1[10,1]));
d1 =(1-abs(d0))*prmtr1[11,1]./(1+abs(prmtr1[11,1]))+ abs(d0)-d0^2;
prm_new[10,1]=2*d0;
prm_new[11,1]=-1*(d0^2 + d1) ;

retp(prm_new);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc V_PROB(EV,HE);  @ CALCULATES Pr[Yt/St,Yt-1]@
local VAL;

VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@

proc (7) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor0,xfor1,xfor2,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,confint,lowerCI,upperCI,stderr,stdprederr,xforecast,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,B_M,PR_TTM,PR_TLM,twostep,
pr_tl00,pr_tl01, pr_tl02, pr_tl10,pr_tl11,pr_tl12,pr_tl20,pr_tl21,pr_tl22,
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22,
pi2hat_00,pi2hat_01,pi2hat_02,pi2hat_10,pi2hat_11,pi2hat_12,
pi2hat_20,pi2hat_21,pi2hat_22,B_T2L00,B_T2L01,B_T2L02,
B_T2L10,B_T2L11,B_T2L12,B_T2L20,B_T2L21,B_T2L22;  

EXTERNAL PROC TRANS, V_PROB;
    
    PRMTR=TRANS(PRMTR1);
 
    P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
    P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
    P22=1;             @Pr[St=2/St-1=2]@
                
    siguGM0 = prmtr[3,1];
    siguGM1 = prmtr[4,1];
    sigv= prmtr[5,1]; 
    sigep = prmtr[6:8,1]; 
    rho = prmtr[9,1];
                    
    phi1 = prmtr[10,1];
    phi2 = prmtr[11,1];
    psi=prmtr[12:14,1];    
    k=prmtr[15:17,1];
    
    A=1~0;      @selection matrix@
    H = 1~1~0~0|0~0~1~0;
               
    B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
    P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
    B_LL1=inf[t1]|0|xx[1]|0; 
    P_LL1 = eye(4)*10000;       
    B_LL2=inf[t2]|0|xx[1]|0; 
    P_LL2 = eye(4)*10000;      
                
    subF = phi1~phi2| 1~0; 
    subQ= siguGM0^2~0|0~0;
    vecp = inv(eye(4)-subF.*.subf)*vec(subq);
    vv = reshape(vecp',2,2);
    p_ll0[3:4,3:4] = vv;
    p_ll1[3:4,3:4] = vv;
    p_ll2[3:4,3:4] = vv;

    p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
    p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
    p_ll2[2,2] = sigep[3,1]^2/(1-psi[3]^2);

    PROB__0 = 1; @Unconditional probability for St=0@
    PROB__1 = 0; @Unconditional probability for St=1@
    PROB__2 = 0; @Unconditional probability for St=2@
 
    B_M=ZEROS(T,tvp_no);  @will store (1,1)element of B_TT@
    PR_TTM=ZEROS(T,3);  @will store Pr[S_t=0|Y_t]@
    PR_TLM=ZEROS(T,3);  @will store Pr[S_t=0|Y_{t-1}]@
    stdprederr=zeros(t,1); @standardized prediction error for diagnostic tests@
    confint=zeros(t,2);  @90% confidence interval bands for trend inflation estimates@ 
    xforecast=zeros(t,1); @stores the infinite sum term in the cyclical component@
    twostep = zeros(t,1);
/*  START ITERATION */
/*===================================================================== */

   LIKV=0.0;
   J_ITER = 1;
   do until J_ITER>T;
        
        yt=inf[j_iter]|xx[j_iter];
        dmt = dmat[j_iter];
        sigu=siguGM0*dmt + siguGM1*(1-dmt);
           
        pr_tl00 = PROB__0*p00;         
        pr_tl01 = PROB__0*(1-p00);
        pr_tl02 = 0;
        pr_tl10 = 0;
        pr_tl11 = PROB__1*p11;
        pr_tl12 = PROB__1*(1-p11);
        pr_tl20 = 0;
        pr_tl21 = 0;
        pr_tl22 = p22 *PROB__2;
              
        PR_TLM[J_ITER,.]= (pr_tl00 +pr_tl10+pr_tl20)
                         ~(pr_tl01 +pr_tl11+pr_tl21)
                         ~(pr_tl02 +pr_tl12+pr_tl22); 
 
        Q0=    sigv^2~0~0~0|
               0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
               0~rho*sigu*sigep[1,1]~sigu^2~0|
               0~0~0~0;
        Q1=    sigv^2~0~0~0|
               0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
               0~rho*sigu*sigep[2,1]~sigu^2~0|
               0~0~0~0;
        Q2=    sigv^2~0~0~0|
               0~sigep[3]^2~rho*sigu*sigep[3,1]~0|
               0~rho*sigu*sigep[3,1]~sigu^2~0|
               0~0~0~0;
        F0 =   1~0~0~0 |
               0~psi[1]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F1 =   1~0~0~0 |
               0~psi[2]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F2 =   1~0~0~0 |
               0~psi[3]~0~0|
               0~0~phi1~phi2|
               0~0~1~0;

        FF= phi1~phi2|1~0;
        xfor0 = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
        ome0 =0~0~xfor0| 0~0~0~0;
        xfor1 = k[2] *A*(inv(eye(2)-FF)); 
        ome1 =0~0~xfor1| 0~0~0~0;
        xfor2 = k[3] *A*(inv(eye(2)-FF)); 
        ome2 =0~0~xfor2| 0~0~0~0;

        /* KALMAN FILTER */
        @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
        B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
        B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
        B_TL02 = F2 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=2@
        B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
        B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
        B_TL12 = F2 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=2@
        B_TL20 = F0 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=0@
        B_TL21 = F1 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=1@
        B_TL22 = F2 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=2@
        
        @ P_tl(i,j),  S(t-1)=i, S(t)=j@
        P_TL00 = F0 * P_LL0 * F0' + Q0;
        P_TL01 = F1 * P_LL0 * F1' + Q1;
        P_TL02 = F2 * P_LL0 * F2' + Q2;
        P_TL10 = F0 * P_LL1 * F0' + Q0;
        P_TL11 = F1 * P_LL1 * F1' + Q1;
        P_TL12 = F2 * P_LL1 * F2' + Q2;
        P_TL20 = F0 * P_LL2 * F0' + Q0;
        P_TL21 = F1 * P_LL2 * F1' + Q1;
        P_TL22 = F2 * P_LL2 * F2' + Q2;
 
        @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
        F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
        F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
        F_CAST02= yt- H * B_TL02 - ome2* B_TL02;
        F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
        F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
        F_CAST12= yt- H * B_TL12 - ome2* B_TL12;
        F_CAST20= yt- H * B_TL20 - ome0* B_TL20;
        F_CAST21= yt- H * B_TL21 - ome1* B_TL21;
        F_CAST22= yt- H * B_TL22 - ome2* B_TL22;
 
        @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
        SS00= H * P_TL00 * H' ; 
        SS01= H*  P_TL01 * H' ;
        SS02= H*  P_TL02 * H' ;
        SS10= H * P_TL10 * H' ;
        SS11= H * P_TL11 * H' ;
        SS12= H * P_TL12 * H' ;
        SS20= H * P_TL20 * H' ;
        SS21= H * P_TL21 * H' ;
        SS22= H * P_TL22 * H' ;
 
        /*Updating*/
        B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
        B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
        B_TT02 = B_TL02 + (P_TL02 * H' *inv(SS02)) * F_CAST02;
        B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
        B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
        B_TT12 = B_TL12 + (P_TL12 * H' *inv(SS12)) * F_CAST12;
        B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
        B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
        B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

        P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
        P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
        P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
        P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
        P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
        P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
        P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
        P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
        P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;
           
       /* HAMILTON FILTER */   

        PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
        PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
        PR_VL02=0;
        PR_VL10=0;
        PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
        PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
        PR_VL20=0;
        PR_VL21=0;
        PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
        
        @CONDITIONAL DENSITY TIMES WEIGHT@

        PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                                  
        @Pr[St,St-1/Yt]@
        PRO_00=PR_VL00/PR_VAL;       
        PRO_01=PR_VL01/PR_VAL;
        PRO_02=PR_VL02/PR_VAL;
        PRO_10=PR_VL10/PR_VAL;
        PRO_11=PR_VL11/PR_VAL;
        PRO_12=PR_VL12/PR_VAL;
        PRO_20=PR_VL20/PR_VAL;
        PRO_21=PR_VL21/PR_VAL;
        PRO_22=PR_VL22/PR_VAL;
        
        PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
        PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
        PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@
        PR_TTM[J_ITER,.]=PROB__0~prob__1~prob__2;
       /* COLLAPSING TERMS */
      
       if(prob__0!=0);
       B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
       endif;
       if(prob__1!=0);
       B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
       endif;
       if(prob__2!=0);
       B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
       endif;
       if(prob__0!=0);
       P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
       PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
       PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
       endif;
       if(prob__1!=0);
       P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
       PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
       PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
       endif;
       if(prob__2!=0);
       P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
       PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
       PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
       endif;
 
       @Estimates of the unobserved components@
       B_TT=   B_TT00*PRO_00+B_TT01*PRO_01+B_TT02*PRO_02+
                B_TT10*PRO_10+B_TT11*PRO_11+B_TT12*PRO_12+
                B_TT20*PRO_20+B_TT21*PRO_21+B_TT22*PRO_22;
       B_M[J_ITER,.]=B_TT';

       B_T2L00 =  F0*F0 *B_LL0;  
       B_T2L01 = F1*F1 *B_LL0;
       B_T2L02 = F2*F2* B_LL0;
       B_T2L10 = F0*F0* B_LL1;
       B_T2L11 = F1*F1* B_LL1;
       B_T2L12 = F2*F2* B_LL1;
       B_T2L20 =F0*F0* B_LL2;
       B_T2L21 = F1*F1* B_LL2;
       B_T2L22 =  F2*F2*B_LL2;
       
       /*Calculating predicted one quarter ahead inflation from the UC model*/
       pi2hat_00 = H * B_T2L00 + ome0* B_T2L00;
       pi2hat_01= H * B_T2L01 + ome1* B_T2L01;
       pi2hat_02 = H * B_T2L02 + ome2* B_T2L02;
       pi2hat_10 = H * B_T2L10 + ome0* B_T2L10;
       pi2hat_11 = H * B_T2L11 + ome1* B_T2L11;
       pi2hat_12 = H * B_T2L12 + ome2* B_T2L12;
       pi2hat_20 = H * B_T2L20 + ome0* B_T2L20;
       pi2hat_21 = H * B_T2L21 + ome1* B_T2L21;
       pi2hat_22 = H * B_T2L22 + ome2* B_T2L22;

       twostep[j_iter,.]= (pi2hat_00[1,1])*pr_tl00 +  pi2hat_01[1,1]*pr_tl01+ 
                pi2hat_02[1,1]*pr_tl02+ pi2hat_10[1,1]*pr_tl10+
                pi2hat_11[1,1]*pr_tl11+ pi2hat_12[1,1]*pr_tl12+ 
                pi2hat_20[1,1]*pr_tl20+ pi2hat_21[1,1]*pr_tl21+
                pi2hat_22[1,1]*pr_tl22;
 
       /*Calculating 90% Confidence bands for trend inflation estimates*/
       lowerCI=(B_TT00[1,1]-1.645*sqrt(P_TT00[1,1]))*PR_tl00+
               (B_TT01[1,1]-1.645*sqrt(P_TT01[1,1]))*PR_tl01+
               (B_TT02[1,1]-1.645*sqrt(P_TT02[1,1]))*PR_tl02+
               (B_TT10[1,1]-1.645*sqrt(P_TT10[1,1]))*PR_tl10+
               (B_TT11[1,1]-1.645*sqrt(P_TT11[1,1]))*PR_tl11+
               (B_TT12[1,1]-1.645*sqrt(P_TT12[1,1]))*PR_tl12+
               (B_TT20[1,1]-1.645*sqrt(P_TT20[1,1]))*PR_tl20+
               (B_TT21[1,1]-1.645*sqrt(P_TT21[1,1]))*PR_tl21+
               (B_TT22[1,1]-1.645*sqrt(P_TT22[1,1]))*PR_tl22;
       upperCI=(B_TT00[1,1]+1.645*sqrt(P_TT00[1,1]))*PR_tl00+
               (B_TT01[1,1]+1.645*sqrt(P_TT01[1,1]))*PR_tl01+
               (B_TT02[1,1]+1.645*sqrt(P_TT02[1,1]))*PR_tl02+
               (B_TT10[1,1]+1.645*sqrt(P_TT10[1,1]))*PR_tl10+
               (B_TT11[1,1]+1.645*sqrt(P_TT11[1,1]))*PR_tl11+
               (B_TT12[1,1]+1.645*sqrt(P_TT12[1,1]))*PR_tl12+
               (B_TT20[1,1]+1.645*sqrt(P_TT20[1,1]))*PR_tl20+
               (B_TT21[1,1]+1.645*sqrt(P_TT21[1,1]))*PR_tl21+
               (B_TT22[1,1]+1.645*sqrt(P_TT22[1,1]))*PR_tl22;
       confint[j_iter,.] = lowerCI~upperCI;

       /*Calculating standardized prediction errors for diagnostic tests*/
       stderr= (f_cast00[1,1]/sqrt(ss00[1,1]))*pr_tl00+ (f_cast01[1,1]/sqrt(ss01[1,1]))*pr_tl01+
               (f_cast02[1,1]/sqrt(ss02[1,1]))*pr_tl02+ (f_cast10[1,1]/sqrt(ss10[1,1]))*pr_tl10+
               (f_cast11[1,1]/sqrt(ss11[1,1]))*pr_tl11+ (f_cast12[1,1]/sqrt(ss12[1,1]))*pr_tl12+ 
               (f_cast20[1,1]/sqrt(ss20[1,1]))*pr_tl20+ (f_cast21[1,1]/sqrt(ss21[1,1]))*pr_tl21+
               (f_cast22[1,1]/sqrt(ss22[1,1]))*pr_tl22;
       STDPREDERR[j_iter,.] = stderr;
       
       /*Calculating the infinite sum term in the cyclical component for inflation*/
       xforecast[j_iter,.]= (xfor0 *B_TL00[3:4,1])*pr_tl00 +
                            (xfor1 *B_TL01[3:4,1])*pr_tl01 +
                            (xfor2 *B_TL02[3:4,1])*pr_tl02 +
                            (xfor0 *B_TL10[3:4,1])*pr_tl10 +
                            (xfor1 *B_TL11[3:4,1])*pr_tl11 +
                            (xfor2 *B_TL12[3:4,1])*pr_tl12 + 
                            (xfor0 *B_TL20[3:4,1])*pr_tl20 +
                            (xfor1 *B_TL21[3:4,1])*pr_tl21 +
                            (xfor2 *B_TL22[3:4,1])*pr_tl22;

       if J_ITER < START; goto skip; endif;
       LIKV = LIKV -LN(PR_VAL);
       skip:
       J_ITER = J_ITER+1;
       endo;
  
/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],PR_TTM[START:T,.],PR_TLM[START:T,.],stdprederr[start:t,.],xforecast[start:t,.],confint[start:t,.],twostep[start:t,.]);
ENDP;


@==========================================================@

@Kim's smoothing@
proc SMOOTH(prtt,prtl,prmtr0);

local   pr_sm1,pr_sm2,pr_sm3, j_iter, TT,prmtr, pr_tr, p11,p12,p13,p21,p22,p23,p31,p32,p33,
        pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3,pr_sm11,pr_sm12,pr_sm13,pr_sm21,pr_sm22,pr_sm23,pr_sm31,pr_sm32,pr_sm33;

      TT=ROWS(PRTT);
      prmtr=prmtr0;
      p11=prmtr[1,1]; 
      p12=1-p11; 
      p13=0; 
      p21=0; 
      p22=prmtr[2,1]; 
      p23=1-p22;
      p31=0; 
      p32=0; 
      p33=1;
 
      pr_tt1=prtt[.,1];
      pr_tt2=prtt[.,2];
      pr_tt3=prtt[.,3];
      pr_tl1=prtl[.,1];
      pr_tl2=prtl[.,2];
      pr_tl3=prtl[.,3];
           
      pr_sm1=pr_tt1;      @ pr_sm1 will contain Pr[S_t=1|Y_T]@
      pr_sm2=pr_tt2;      @ pr_sm2 will contain Pr[S_t=2|Y_T]@
      pr_sm3=pr_tt3;      @ pr_sm3 will contain Pr[S_t=3|Y_T]@
          

     j_iter=TT-1;
     do until j_iter < 1;
  
     if pr_tl1[j_iter+1,1]<0.000001; pr_tl1[j_iter+1,1]=0.0000001; endif;
     if pr_tl2[j_iter+1,1]<0.000001; pr_tl2[j_iter+1,1]=0.0000001; endif;
     if pr_tl3[j_iter+1,1]<0.000001; pr_tl3[j_iter+1,1]=0.0000001; endif;
 
        pr_sm11=pr_sm1[j_iter+1,1]*p11*  pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm12=pr_sm2[j_iter+1,1]*p12*  pr_tt1[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm13=pr_sm3[j_iter+1,1]*p13*  pr_tt1[j_iter,1]/ pr_tl3[j_iter+1,1];
 
        pr_sm21=pr_sm1[j_iter+1,1]*p21*  pr_tt2[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm22=pr_sm2[j_iter+1,1]*p22*  pr_tt2[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm23=pr_sm3[j_iter+1,1]*p23*  pr_tt2[j_iter,1]/ pr_tl3[j_iter+1,1];
      
        pr_sm31=pr_sm1[j_iter+1,1]*p31*  pr_tt3[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm32=pr_sm2[j_iter+1,1]*p32*  pr_tt3[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm33=pr_sm3[j_iter+1,1]*p33*  pr_tt3[j_iter,1]/ pr_tl3[j_iter+1,1];
       
        pr_sm1[j_iter,1]=pr_sm11+pr_sm12+pr_sm13;
        pr_sm2[j_iter,1]=pr_sm21+pr_sm22+pr_sm23;
        pr_sm3[j_iter,1]=pr_sm31+pr_sm32+pr_sm33;
    
     j_iter=j_iter -1;
     endo;


retp(pr_sm1~pr_sm2~pr_sm3);  
endp;

/*=================================================================================================================*/
 